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By combining the Baeriswyl wavefunction with equilibrium and time-dependent variational prin¬ 
ciples, we develop a non-equilibrium formalism to study quantum quenches for two dimensional spin¬ 
less fermions with nearest-neighbour hopping and repulsion. The variational ground state energy 
and the short time dynamics agree convincingly with the results of numerically exact simulations. 
We find that depending on the initial and final interaction strength, the quenched system either 
exhibits undamped oscillations or relaxes to a time independent steady state. The time averaged 
expectation value of the CDW order parameter rises sharply when crossing from the steady state 
regime to the oscillating regime, indicating that the system, being non-integrable, shows signs of 
thermalization with an effective temperature above or below the equilibrium critical temperature, 
respectively. 

PACS numbers: 71.10.Fd,05.30.Fk,05.70.Ln 


Introduction — The past decade has witnessed a 
massive resurgence of interest in the coherent dynam¬ 
ics of quantum many-body systems far from equilibrium 
[H-Q. This activity has been catalyzed mostly by de¬ 
velopments in cold atom experiments ii , which allow 
for long coherence and explicitly tracking the evolution 
of well-isolated many-body systems in real time. At the 
same time, there is also a revival of interest in the ultra¬ 
fast time evolution of solid-state systems. A ubiquitous 
paradigm in the theoretical study of non-equilibrium dy¬ 
namics is the so-called quantum quench, where one starts 
from the ground state of a Hamiltonian, and then in¬ 
stantaneously changes a parameter, so that the system 
wavefunction evolves under a different Hamiltonian. 


In contrast to the one-dimensional (ID) case, there are 
only few exact theoretical methods available to treat real¬ 
time evolution in 2D quantum systems. For example, the 
DMRG technique [5-ll| is challenging in 2D, there are 
almost no 2D situations that are integrable via the Bethe 
ansatz, and exact numerical diagonalization is faced with 
a severely rapid growth of Hilbert space size. There are 
thus relatively few studies of interaction quenches in 2D, 
e.g.. Refs, [ll, 1^ " 
namics, or Ref. [it 
interactions. 


focusing mainly on short time dy- 
on time evolution in the absence of 


Variational wavefunctions are vital in equilibrium con¬ 
densed matter physics, e.g., the variational BCS wave- 
function for superconductivity, the Yosida wavefunction 
for the Kondo model 1^, and the Gutzwiller wavefunc¬ 
tion (GWF) [3i [l3| for the Hubbard model. The pa¬ 
rameters of a well-chosen variational wavefunction can 
incorporate the most essential physics of a complex 
many-body system. In contrast to equilibrium, the non¬ 
equilibrium variational principle [l^ has been used less 
extensively in correlated systems. The bosonic version of 


the GWF was often used during the last decade to ex¬ 
plore dynamics in the Bose-Hubbard model i 3“21|. Re¬ 
cently, computational methods using Monte Carlo eval¬ 
uation were formulated for the evolution of variational 
wavefunctions with large numbers of parameters 22 - 2^ . 
In additiom the GWF within the Gutzwiller approxima¬ 
tion H, [iTj was recently used to study the time evolu¬ 
tion 23 of the fermionic Hubbard model. This approx¬ 
imation gives the exact solution of the GWF in infinite 
dimensions 12611. 


The GWF [l6,[l7| consists of a non-interacting wave- 
function acted on by an operator which projects out 
states not favoured by the interaction (i.e. double oc¬ 
cupation for the spinful Hubbard model). A completely 
opposite variational approach, the Baeriswyl wavefunc¬ 
tion (BWF) [ 23 , dll, is based on a fully localized wave- 
function (an eigenfunction of the interaction term) and 
is acted on by a kinetic energy projector. The projector 
has the effect of promoting the hopping of particles in an 
otherwise completely localized system. In both of these 
approaches, the strength (coefficient) of the projection 
operator serves as the variational parameter. 


In the present work, we will use the BWF to describe a 
model of 2D spinless fermions on the square lattice with 
nearest neighbor repulsion at half filling. The reference 
state is a charge density wave (GDW) with one of the 
sublattices exactly filled. An advantage of this approach, 
as we show, is that it is possible to perform the calcu¬ 
lations exactly within the variational manifold, without 
the requirement of additional approximations. In addi¬ 
tion, the wavefunction is exact at the two extreme limits 
of infinite and zero interaction. We show that equilibrium 
properties are well-reproduced at all interactions. We use 
the time-dependent variational principle [3,0| to inves¬ 
tigate the dynamics after interaction quenches, and show 
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that the system either relaxes to a steady state or os¬ 
cillates indefinitely. We map out the “quantum quench 
phase diagram” according to this criterion and conjecture 
an interpretation based on the thermalization tempera¬ 
ture. The analysis involves characteristics of the CDW 
order parameter. Density waves have long been a cen¬ 
tral topic in the study of electronic phases in the solid 
state [30j. Recent experiments have also addressed their 
non-equilibrium dynamics for CDW jsij as well as of re¬ 
lated spin density (antiferromagnetic) patterns 32, 

The present work studies CDW dynamics for a 2D inter¬ 
acting system, a topic which at present lies outside the 
reach of exact methods. 

Model — The spinless fermionic Hubbard model is 
defined on the 2D square lattice as H = Hkin + Hint 
with 


Hkin — 2 'y y 

{n,m) 


TT _Y_ 

^int — 2 / ^ 

{n,m) 


( 1 ) 


where = c^Cn is the particle density operator at site 
n, and c„ annihilates a fermion from site n. The sum 
over (n, m) runs through all lattice sites and keeps only 
nearest-neighbor pairs. J is the strength of the single¬ 
particle hopping integral, used as the unit of energy in 
the following, thus it is suppressed, H > 0 is the strength 
of the nearest neighbor Coulomb repulsion. The BWF 
amounts to starting from the H > oo state and using 
the kinetic energy as a projection operator 

I'f's) = Ns^e^p{aHk,n) |CDW), (2) 


where |CDW) is the charge ordered ground state H in 
the atomic limit, where the charges follow a checkerboard 
pattern on the lattice sites (see Fig. [T]), Nb is the overall 
normalization factor. 

The BWF, due to its CDW ground state [s^ pos¬ 
sesses an appealing wavefunction in momentum space as 
well. The fully polarized CDW wavefunction is |CDW) = 

n c^-q) |0)^ ’^'^ere Q = (7r/a,7r/a) is the 

kG R.BZ 

ordering wavevector, a the lattice constant, the prod¬ 
uct goes through the reduced (magnetic) Brillouin zone 
(RBZ) and |0) is the fermionic vacuum. The normalized 
variational wavefunction assumes the form 


l^s) 


n 

kGRBZ 


exp[5e(k)]cj^ -I- exp[-Q!e(k)]c;^__Q 
y/2 cosh[2Rede(k)] 


|0) (3) 


with e(k) = — cos{kxa) — cos{kya). Here, a = a + ir] is 
the complex variational parameter. Its imaginary part 
becomes relevant when studying the quench problem. 

Equilibrium properties — The variational problem 
is solved by minimizing the ground state energy with 
respect to the variational parameter a as Eq = 
min('I'B|i7|'I'B), and our conventions imply that negative 

a 

real values of a minimize the above functional, namely 



FIG. 1. (Color online) The variational ground state energy, 
normalized by its non-interacting value (blue solid line), is 
compared to 2D DMRG data (red circles) and ED on 4 x 4 
cluster (black squares). The lower inset visualizes the CDW 
order parameter uq (red dashed line), the equal time struc¬ 
ture factor from Eq. 0 for momentum Q (green dash-dotted 
line) and the variational parameter, a (blue line) as a function 
of the interaction strength. The upper inset depicts a small 
segment of the fully polarized CDW state. 


a < 0 and rj = 0. The expectation value of the kinetic 
and interaction energy is evaluated exactly with the vari¬ 
ational wavefunction as {Hkin) = e{ 0 )l 2 and 


= -nv'M + 


(4) 


2 =1,2,3 


where N is the total number of lattice sites and [34 1 




cos [2776 (k)] 

2 cosh[2ae(k)] 


T- X—' C0s(/Ca;0) - ,, 

72 =^- 7 .—- tanh[2ae(k)]. 


^3 = E 


cos{kxa) sin[277e(k)] 
2 cosh[2ae(k)] 


(5a) 

(5b) 

(5c) 


Note that we obtain an exact, closed expression for the 
variational energy, {Hkin) + {Hint) hr 2D for this model. 
So far, closed expressions have not been derived based 
on the GWF for any two-dimensional model. Moreover, 
the energy expectation values are valid also for hyper- 
cubic lattices in arbitrary dimension d, after making the 
replacement e(k) = — cos(fcia), and the k sums run 
over d wavevector components. 

The occupation number of momentum states is rik = 
(CkCk) = 1/{1 + exp[—4ae(k)]}. This expression corre¬ 
sponds to the occupation number of a non-interacting 
system at finite temperature with T = —l/4a. In the 
following we will make use of this interpretation. The ini¬ 
tial CDW state thus possesses infinite temperature with 
rik = 1/2, as follows from |CDW) as well, albeit it differs 



























3 



FIG. 2. (Color online) The quantum quench phase dia¬ 
gram for Vi and Vf is shown. For small Vf, a time inde¬ 
pendent steady state is reached which could be interpreted 
as an indicator thermalization at an effective temperature 
higher than the equilibrium CDW transition temperature 
with r]{t oo) ~ t and a{t) saturating to a finite value, 
as depicted in the lower inset for Vi = 1.5, Vf = 0.25 and 
0 < f < 200. For sizable quenches, on the other hand, the 
system periodically returns to its initial state, visualized in 
the upper inset for Vi = 1.5, Vf = 1. For Vi = oo, the transi¬ 
tion from the steady state to oscillating behaviour occurs at 
Vf « 0.97. 


from a trivial T = oo state due to the anomalous cor¬ 
relations with wavevector Q from the CDW. The charge 
density at lattice site R is 

^+«QCOs(QR), ng = ^, (6) 

with UQ the CDW order parameter, describing the Q 
periodic charge oscillations. 

The Qth Fourier component of the equal-time con¬ 
nected structure factor (SF), which measures non-local 
properties and CDW correlations, is 


the ground state energy, the CDW order parameter ng 
and the optimal variational parameter are shown in Fig. 
[TJ together with numerical data using 2D DMRG and 
exact diagonalization (ED) on 4 x 4 cluster. The 2D 
DMRG data is obtained by extrapolating the energies for 
infinitely long cylinders with circumferences L = 6,8,10 
to the thermodynamic limit. The agreement between the 
variational and numerically obtained energies is remark¬ 
able. The ED matches the DMRG data for larger V, 
indicating short correlation length, such that the rela¬ 
tively small system size does not influence the ground 
state energy. 

The system is always in a CDW state, except for 
a —>• — oo, corresponding to the non-interacting limit 
E = 0. This is in agreement with what is expected for 
the spinless Hubbard model on the square lattice 
The CDW phase appears because of the square shaped 
Fermi surface at half Hlling, in the case of perfect nest¬ 
ing [ 3 ^, and also due to the log-divergent density of 
states upon approaching half filling as ^ ln(l/u;). 

Quantum quench — We now turn to the investigation 
of the quantum quench, when an initial E is changed 
suddenly to Vf at t = 0. The time-dependent wave- 
function is of Baeriswyl form, |'I'B(t)) as well, and the 
quench amounts to allowing for time-dependent varia¬ 
tional parameters a{t) and r](t) in Eq. ([3|). Their time- 
dependence follows from the time-dependent variational 
principle, which requires the minimization of the real La- 
grangian, dehned as 


L{t) = -Im({4'B|9t4'B)) - (4 'b|R|4'b), 


( 8 ) 


with respect to the time-dependent variational param¬ 
eters. Using L{t) = (Hkin) {dtV- 1) - (Rmt) with the 
expectation values taken from Eq. with the time- 
dependent variational parameters, we finally arrive to 
the Euler-Lagrange equations, determining the quench 
dynamics as 


5'c(Q) = (^'bIpqP-qI'I'b) - Nlnql'^ = 

1 1 ^ cos^[ 2 ? 7 e(k)] 

2 2N ^ cosh^[2ae(k)] ’ 


dta = 


djHint) 

dr] 

djHkiu) 

da 


dtr] = l + 


djHint) 

da 

djHkin) 

da 


(9) 


where pq = Sk'^k ^k-q, and gives 0 in the fully 
polarized CDW state and 1/2 in the absence of CDW 
correlations, corresponding to infinite and zero effective 
temperatures. Therefore, the value of the SF is related to 
the amount of CDW order and its effective temperature 
of the system. It is time independent in equilibrium, 
shown in Fig. [TJ and reflects the behaviour of ng. 

The variational wavefunction is exact in two opposite 
limits, V = 0 ja = — 00 ) and E = 00 (a = 0). In the 
former case, the ground state energy is Eq = — 4 A^/ 7 r^, 
while for the latter, the kinetic energy is suppressed by 
the CDW and the interaction energy reaches its mini¬ 
mum, therefore Eq = 0. In between these two extrema. 


supplemented with the initial conditions Q!( 0 ) = at, orig¬ 
inating from the equilibrium configuration for E, and 
77 ( 0 ) = 0. From these, the total energy is conserved after 
the quench, as expected. 

From Eq. m, the two limiting cases are recovered. 
First, when no quench was performed, ajt) = at and 
77 (f) = 0. Second, in the case of quenching from E = 00 
to Vf = 0, the at = 0 initial condition fixes ajt) = 0, 
while 77 (f) = f. For general values of initial and final inter¬ 
actions, one needs to i nteg rate these coupled differential 
equations numerically [3^. The obtained phase diagram 
is shown in Fig. |2| Eq. (|T|) in 2D is non-integrable, there¬ 
fore it is expected to thermalize after a sudden quench. 
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FIG. 3. (Color online) Top left: The short time decay of 
the CDW order parameter from time-dependent 2D DMRG 
for an infinite cylinder of circumference L = 8 (solid lines) 
and variational calculation (dashed lines) for Vi = oo and 
Vf = 0 (blue), 1 (red) and 2 (black). Top right: The typi¬ 
cal behaviour of the CDW order parameter is plotted in the 
steady state region (blue line), vanishing as 1/t, for Vi = 1.5 
and Vf = 0.25 and in the oscillatory region (black line) for 
Vi = 1.5 and Vf = 1. Bottom: The time average value of the 
CDW order parameter (blue solid line), the amplitude (red 
dotted line) and the frequency (black dash-dotted line) of os¬ 
cillations, together with the structure factor for momentum 
Q (green dashed line) are shown for Vi = 1.5. 


and we argue that our simple variational scheme is able 
to capture some of its physics. 

The time dependence of the CDW order parame¬ 
ter is calculated from Eq. (|6]) after inserting the time- 
dependent variational parameters into If. The short time 
behaviour from the BWF agrees well with that from 2D 
time-dependent DMRG obtained using the algorithm in¬ 
troduced in Ref. ll|, as shown in Fig. [3] The Vi = oo ^ 
Vf = 0 quench is exact within the variational framework, 
which also allows us to check the temporal validity of the 
numerics. 

For longer times, the variational solution either oscil¬ 
lates around its time average, or reaches a steady state 
after the quench [s^. For quenches with small Vf, a 
time independent steady state is reached with r]{t) in¬ 
creasing linearly with time and a(t) saturating to a fixed 
value. This can be thought of as the heating up of the 
system, such that it eventually thermalizes to a higher 
effective temperature than the equilibrium CDW transi¬ 
tion temperature, thus CDW is absent. The CDW order 
parameter decays to zero since due to rjit) ^ t, the nu¬ 
merator in If oscillates fast and kills the integral with in¬ 
creasing t. In particular, for arbitrary spatial dimension 
d, a riQ^t) ^ decay is found, in accord with special 

quenches of the ID Heisenberg chain [^. A similar decay 
is found for the antiferromagnetic order of the 2D spinful 
Hubbard model [l^ . quenched to a non-interacting sys¬ 


tem. While in ID, the CDW order parameter oscillates 
and changes also sign after the quench IJ, [37| , it only 
oscillates around its envelope function in 2D but does not 
change sign, as shown in Fig. [3l This is supported by 
the exact result when quenching from V) = oo to V/ = 0, 
nQ(t) = Jg(2t)/2, valid in arbitrary dimension d, and 
Jo (a;) is the zeroth Bessel function of the first kind, in 
agreement with our numerical findings. 

For larger Vf, on the other hand, an oscillating solu¬ 
tion is found for both a(t) and rj{t), and the system pe¬ 
riodically returns to its initial state, and seemingly does 
to reach a time independent steady state, illustrated in 
Fig. [21 Within our approach, we interpret this as an 
indicator of thermalization to a thermal state with lower 
effective temperature than the equilibrium CDW tran¬ 
sition temperature, hence the resulting state possesses 
a finite CDW order. The oscillations are artifacts of 
the variational calculations, but the value it oscillates 
around could correspond to that in the thermalized state. 
Throughout the time evolution, a{t) stays negative, and 
its equilibrium interpretation as an effective temperature 
through the momentum distribution remains true after 
the quench, with a time-dependent effective temperature 
T{t) = -l/4a(t). 

In order to further characterize the CDW phase, 
we consider the time average of the order param¬ 
eter as ng = limt_>.oo \ Jq nQ(t')dt', together with 
the amplitude, Ang and frequency, wg of the oscil¬ 
lations of nQ[t) in the oscillating regime, shown in 
Fig. [3] for the representative case oi Vi = 1.5. Here, 
Ang = limt_>oo {inax[ng(t)] — min[ng(t)]} in the oscil¬ 
lating regime and wg is the basic harmonics of the oscil¬ 
lations in ngit), evaluated from Fourier analysis. Due to 
the non-integrability of the model, the system is expected 
to thermalize Q. For small Vf, the CDW gap and transi¬ 
tion temperature is small and the effective temperature, 
the system reaches, is above the equilibrium CDW tran¬ 
sition temperature, therefore the CDW order is absent in 
this case, as evidenced in Fig. [3] For larger final inter¬ 
action, however, the corresponding equilibrium transition 
temperature increases, and becomes equal to the effective 
temperature of the thermalized system, which marks the 
sharp rise of ng. From that point on, with increasing 
Vf, the system’s effective temperature is always smaller 
than the equilibrium transition temperature, therefore 
the system exhibits long range CDW order, mimicking a 
thermalized state. This scenario is further corroborated 
by focusing on other characteristics of the order, such 
as the amplitude and frequency of the oscillations in the 
oscillating regime. 

The equal time Qth structure factor after the quench 
is obtained by inserting the time-dependent variational 
parameters in Eq. ©• Its long time average, 5'c(Q) is 
plotted in Fig. [31 At the critical Vf, where the transition 
occurs, the variational wavefunction is almost metallic, 
and contains enhanced CDW correlations away from this 
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point. However, only strong interaction can profit from 
these correlations and drive the dynamical CDW transi¬ 
tion. The effective temperature after the quench is ex¬ 
pected to increase monotonically with \ Vf — Vi\, while the 
equilibrium CDW gap grows with Vf. When these two 
energy scales become comparable, the dynamical phase 
transition occurs.. 

Conclusion — We have studied quantum quenches 
for 2D spinless fermions on the square lattice with 
nearest-neighbour hopping and repulsion using the vari¬ 
ational BWF, which are unaccessible by exact methods. 
Depending on the initial and final interaction strength, 
the system reaches either a time independent steady state 
or oscillates periodically in time. We argue by investigat¬ 
ing the characteristics of the CDW order parameter and 
the equal time structure factor that the former and latter 
behave similarly to what is expected from a thermal state 
with effective temperature larger and smaller than the 
equilibrium CDW transition temperature, respectively. 

Our work opens up a number of interesting questions 
worth pursuing in the future, such as considering bosons 
instead of fermions, the effect of disorder, the influence 
of imperfect nesting, incorporating the spin degree of 
freedom, as well as the improvement of the variational 
method by introducing variational parameters for each k 
mode. 
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SUPPLEMENTARY MATERIAL FOR 
’’QUANTUM QUENCH IN 2D USING THE 
VARIATIONAL BAERISWYL WAVEFUNCTION” 
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momentum space as 

Hkin = ^ e(k)cjck, (SI) 

k 

(S2) 

kjk'.q 

The interaction term can be decoupled using Wick’s theo¬ 
rem, yielding Eq. O- The first term in Eq. ® is the con¬ 
ventional Hartree term, Ii also comes from the Hartree 
decoupling by taking the anomalous, (cj^Ck-q) 7 ^ 0 ex¬ 
pectation values into account. The Fock terms give rise 
to / 2 , 3 , containing normal, (cj^Ck) and anomalous expec¬ 
tation values, respectively. Due to the rotational sym¬ 
metry of the square lattice, the x and y directions are 
completely equivalent to each other, hence the cos{kxa) 
factor in 12 ^ 3 , coming from the q dependence of the in¬ 
teraction, together with the appropriate combinatorial 
factors, gives the total variational energy. While the in¬ 
teraction term is even in the variational parameters, the 
kinetic energy is odd in a and independent of rj, therefore 
their balance determines the optimal variational param¬ 
eter. 


TIME DEPENDENCE OF THE VARIATIONAL 
PARAMETERS 


FIG. SI. (Color online) The time dependence of the vari¬ 
ational parameters is plotted for Vi = 1.5 and several final 
Vf. The critical interaction strength, separating the steady 
state and oscillating regimes is around Vf ~ 0.591. Note the 
different horizontal and vertical scales. 


ENERGY EXPECTATION VALUE 

Due to the structure of the wavefunction, the problem 
is more tractable when writing the Hamiltonian H in 


The numerical solution of Eqs. yield the time- 
dependent variational parameters, whose behaviour is 
plotted in Fig. [STl In the steady state regime, the ri{t) 
keeps on increasing linearly with time, while a{t) reaches 
a time independent steady state value. However, the re¬ 
laxation time required to reach this increases with in¬ 
creasing Vf. At a critical final interaction strength, the 
solution changes abruptly and both parameters exhibit 
periodic oscillations. The oscillation frequency increases 
sharply with V/, and contain higher harmonics as well 
close to the critical interaction strength. With increas¬ 
ing Vf, however, a single frequency describes more and 
more reliably the data. 

















